1. Field of the Invention
This invention relates to radiation therapy. Embodiments of the invention are implemented as a software system, called NXEGS, that uses Monte Carlo simulation of radiation transport for radiotherapy treatment planning (RTP) and includes the following capabilities: commissioning of a linear accelerator to construct an equivalent source model, simulation of radiation transport through beam modifiers, and simulation of radiation transport in a patient or phantom to obtain the resulting dose distribution. Embodiments of the invention are implemented using a digital computer programmed to execute the NXEGS software and algorithms described herein.
2. Description of Related Art
The following background citations will be referred to by the references indicated in brackets. The entirety of all of these citations are incorporated herein by reference.
[1] W. R. Nelson and D. W. O. Rogers, xe2x80x9cStructure and Operation of the EGS4 Code Systemxe2x80x9d in Monte Carlo Transport of Electrons and Photons, T. M. Jenkins, W. R. Nelson and A. Rindi, eds., Plenum, 287-305, (1989).
[2] D. W. O. Rogers and A. F. Bielajew, xe2x80x9cCode Accuracyxe2x80x9d (Section III, pp. 492-522) in xe2x80x9cMonte Carlo techniques of electron and photon transport for radiation dosimetryxe2x80x9d in The Dosimetry of Ionizing Radiation, Vol. III, K. Kase, B. Bjangard and F. Attix, eds., Academic Press, 427-539, (1990).
[3] D. W. O. Rogers, B. A Faddegon, G. X. Ding, C. -M. Ma, J. We and T. R. Mackie, xe2x80x9cBEAM: A Monte Carlo code to simulate radiotherapy treatment unitsxe2x80x9d Medical Physics 22 (1995) 503-524.
[4] C. -M. Ma, B. A Faddegon, D. W. O. Rogers and T. R. Mackie, xe2x80x9cAccurate characterization of Monte Carlo calculated electron beams for radiotherapyxe2x80x9d Medical Physics 24 (1997) 401-416.
[5] J. J. DeMarco, T. D. Solberg and J. B. Smathers, xe2x80x9cA CT-based Monte Carlo simulation tool for dosimetry planning and analysisxe2x80x9d Medical Physics 25 (1998) 1-11.
[6] C. L. Hartmann Siantar, et al. xe2x80x9cLawrence Livermore National Laboratory""s PEREGRINE Projectxe2x80x9d UCRL-JC-126732 and in 12th International Conference on the Use of Computers in Radiation Therapy. (1997).
[7] A. F. Bielajew and D. W. O. Rogers, xe2x80x9cVariance-Reduction Techniquesxe2x80x9d in Monte Carlo Transport of Electrons and Photons, T. M. Jenkins, W. R. Nelson and A. Rindi, eds., Plenum, 407-419, 1989.
[8] D. L. Donoho xe2x80x9cDe-noising by soft-thresholdingxe2x80x9d, IEEE Transaction. Information Theory 41 (1995) 613-627.
[9] L. Rudin and S. J. Osher xe2x80x9cTotal variation based image restoration with free local constraintsxe2x80x9d Proceedings of 1st International Conference on Image Processing, Austin, Tex., IEEE Comput. Soc. Press, (1994), p.31-35.
[10] A. L. Ames, D. R. Nadeau and J. L. Moreland (1997) VRML 2.0 Sourcebook Wiley pp. 112-115.
[11] H. A. Bethe (1953) xe2x80x9cMolixc3xa8re""s theory of multiple scatteringxe2x80x9d Phys. Rev. 89, 1256-1266.
[12] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery. xe2x80x9cLinear Regularization Methodsxe2x80x9d in Numerical Recipes in C. The Art of Scientific Computing. Cambridge Press (1992) pp. 309-315, 620-623, 808-815.
[13] xe2x80x9cEvaluated Teletherapy Source Libraryxe2x80x9d (2000) U.S. Pat. No. 6,029,079.
[14] xe2x80x9cCalculation of Radiation Therapy Dose Using All Particle Monte Carlo Transportxe2x80x9d (1999) U.S. Pat. No. 5,870,697.
[15] xe2x80x9cTreatment Planning Method and Apparatus for Radiation Therapyxe2x80x9d (1999) International Patent Application #WO99/40523
[16] J. O. Deasy. xe2x80x9cDenoising of electron beam Monte Carlo dose distributions using digital filtering techniquesxe2x80x9d Phys. Med. Biol. 45 (2000) 1765-1779.
[17] I. Kawrakow and M. Fippel (2000) xe2x80x9cInvestigation of variance reduction techniques for Monte Carlo photon dose calculation using XVMCxe2x80x9d Phys. Med. Biol. 45, 2163-2183.
[18] J. Sempau, S. J. Wilderman and A. F. Bielajew (2000) xe2x80x9cDPM, a fast, accurate Monte Carlo code optimized for photon and electron radiotherapy treatment planning dose calculationsxe2x80x9d Phys. Med. Biol. 45, 2263-2291.
[19] T. Holmes and T. R. Mackie (1994) xe2x80x9cA filtered backprojection dose calculation method for inverse treatment planningxe2x80x9d Med. Phys. 21, 303-313.
[20] J. R. Cunningham (1972) xe2x80x9cScatter-air ratiosxe2x80x9d Phys. Med. Biol. 7, 42-51.
[21] T. R. Mackie, J. W. Scrimger and J. J. Battista (1985) xe2x80x9cA convolution method of calculating dose for 15 MV x-raysxe2x80x9d Med. Phys. 12, 188-196.
[22] M. Kalos and P. Whitlock (1986) Monte Carlo Methods Vol. 1. Wiley. pp. 92-103 and 107-109.
[23] E. Veach and L. J. Guibas (1995) xe2x80x9cOptimally combining sampling techniques for Monte Carlo renderingxe2x80x9d in Computer Graphics Proceedings. SIGGRAPH 95, R. Cook, ed., ACM, 1995. p.419-28.
[24] K. R. Hogstrom, M. D. Mills ad P. R. Almond (1981) xe2x80x9cElectron beam dose calculationsxe2x80x9d Phys. Med. Biol. 26, 445.
[25] B. A. Faddegon and I. Blevis (2000) xe2x80x9cElectron spectra derived from depth dose distributionsxe2x80x9d Med. Phys. 27, 514-526.
[26] A. F. Bielajew and D. W. O. Rogers (1987) xe2x80x9cPRESTA: The Parameter Reduced Electron-Step Transport Algorithm for Electron Monte Carlo Transportxe2x80x9d Nucl. Instr. Meth. Phys. Res. B18, 165-181.
About 1,400,000 people in the US are diagnosed with cancerous tumors each year. Every American has roughly a 50% chance of contracting some form of cancer during their lifetime, and the number of cases is anticipated to rise substantially in the years to come. Roughly 50% of those diagnosed are treated with radiation in the form of high energy electron beams or photon (X-ray) beams.
FIG. 1.1 shows a schematic drawing of the treatment setup. The radiation source is typically a linear accelerator (linac). After it is emitted from the linac, the radiation beam passes through a set of beam modifiers and then on to the patient or phantom. The treatment head, consisting of linac and beam modifiers, is contained in a moveable gantry, whose position and direction are variable. FIG. 1.2 shows more detail of the components within the treatment head. Following the linac, one finds four moveable jaws (xe2x80x9cleftxe2x80x9d and xe2x80x9crightxe2x80x9d, xe2x80x9cupxe2x80x9d and xe2x80x9cdownxe2x80x9d) that determine the width of the photon beam in the plane transverse to the central axis and one or more scrapers used to define the size of the electron beam. Then there are beam modifiers (such as the wedge and port shown in the FIG. 1.2) that are used to modify the intensity profile of the beam.
In order to make the treatment more effective, one must plan quantitatively the exposure to radiation in terms of dose amount, time duration and scheduling of treatments. This requires an accurate simulation of the way in which the radiation beams pass through the patient""s tissue and interact with the tumor. Typically, multiple radiation beams are used on a patient, and their direction and intensity must be chosen carefully to apply adequate dose to the cancer without harming healthy tissues.
Intensity modulated radiation therapy (IMRT) and other new treatment modalities allow much more aggressive treatment therapies. These might prescribe a high dose to a target that is close to a critical tissue, so that relatively small errors in the dose calculation could be of great clinical significance. Monte Carlo simulation is the most accurate method for computation of radiation dose delivery, and so it is desired for radiotherapy treatment planning (RTP).
The Monte Carlo method for radiation therapy involves direct computational simulation of the physics of particle transport in the transport media (the beam modifiers and the patient). Particles interact with the transport media through atomistic processes, the outcomes of which are chosen randomly using scattering cross-sections. xe2x80x9cMonte Carloxe2x80x9d refers to these random choices, which are made with the help of a computerized random number generator.
EGS4
The Monte Carlo computer code EGS4, developed at the Stanford Linear Accelerator (SLAC) and the National Research Council of Canada (NRCC) [1], simulates the radiation transport and deposition of dose within a general geometry, such as a patient or phantom. It has been carefully validated and shown to provide highly accurate dose distributions (e.g. see [2]).
The patient or phantom is divided into spatial cells, or voxels, using a rectangular or nonrectangular grid. A typical cell size is 5 mm in each direction. Within each cell the material properties of the patient are input from CT scan data.
Radiation transport in EGS4 includes three kinds of particles: photons (treated as particles as opposed to waves), electrons and positrons. Particles originate at a source and then travel through the phantom (or patient), interacting with the phantom (or patient) material. These interactions result in scattering of the primary particles and generation of new, secondary particles. Each interaction results in a change in the direction and energy of the scattered particle. The change in energy is due to energy imparted to the new particles and energy deposited in the material.
Each scattering process is described by a scattering cross section, which provides the probability or rate for scattering at a given angle. The relevant scattering processes in the Monte Carlo code, EGS4, are described below. As a particle undergoes scattering, the choice of angle and energy from the scattering cross section is a random choice, which is simulated using a computer generated random number (more correctly called a pseudo-random number).
For each source particle, the resulting trajectories of the primary and the resulting secondary particles is described in EGS4 as a shower. Each shower consists of a primary source particle (either a photon or electron) and all of the secondary particles (including particles generated from secondary particles) that are generated.
A photon moves through a material along a path consisting of straight line segments separated by discrete interactions. In EGS4, four interactions of photons with matter are included: pair production, Compton scattering, the photoelectric effect and Rayleigh scattering. In addition, new photons are produced by interactions of electrons and positrons with matter, including Bremstrahlung emission and annihilation.
The transport of electrons and positrons is more complicated than the transport of photons, because their interaction rate (i.e. cross-section) is much larger. EGS4 is a so-called xe2x80x9cClass IIxe2x80x9d algorithm in which the energy loss and deflection of electrons is split into two components: discrete interactions and continuous energy loss. The continuous energy loss simulates the interactions that involve change of energy below an energy threshold. Interactions involving change of energy above this energy threshold are represented as discrete interactions. These latter include Bremstrahlung emission, Mxc3x6ller scattering, Bhabha scattering and annihilation, and elastic scattering. Along each segment of the path of an electron, the occurrence of one of these processes is determined by sampling from the appropriate cross-section. If a discrete interaction does occur, then the resulting deflection, energy loss and secondary particles are randomly chosen from appropriate cross-sections as well. Energy that is not transferred to the secondary particles is absorbed by the material. When a particle energy falls below a certain cutoff value (ecut and pcut, for electrons and photons respectively), the particle is discarded and all of its energy is deposited into the medium.
Of the electron interactions listed above, elastic scattering is treated differently from the others. It is represented by a multiple-scattering process, which groups many small elastic scatterings into a single multiple-scattering step, as described by Molixc3xa8re""s theory [11]. This is applied at discrete increments along an electron path, with a change of angle in each increment. The transport step is modified to allow for larger steps and for steps that cross a material boundary using the PRESTA method [26]. Since each multiple-scattering interaction corresponds to many elastic scatterings, the energy deposition from multi-scattering is applied continuously along the electron path, according to the xe2x80x9cContinuous Slowing Down Approximationxe2x80x9d (CSDA). This is distinct from the continuous energy deposition used to represent interactions involving small energy change, as described in the previous paragraph.
Once the simulation has finished, then the accumulated energy deposition in each cell of the patient geometry is converted into a dose, which is the ratio of the deposited energy in the cell divided by the total mass in the cell. The resulting dose values have units of energy deposition per mass. The commonly used unit of dose is the gray (denoted Gy), which equals one Joule per kilogram, or the rad (also called a centi-gray or cGy) which is equal to one one-hundreth of a gray.
Accuracy of the computed dose distribution and time required for the computation are two crucial issues for any Monte Carlo method. In general terms, the accuracy level xcex5 and the computational time Tcomp are given by
xcex5=c1Nxe2x88x921/2
Tcomp=c2N=c2c12xcex5xe2x88x922xe2x80x83xe2x80x83(1)
in which N is the number of particles, the constant c1 is a measure of (the square root of) the variance and the constant c2 is a measure of the computational time per history. This shows that the accuracy level and the computational time are closely related.
EGS4 includes a number of variance reduction methods that reduce the error xcex5 for a given number N of particles by reducing the variance c12 [7]. A variance reduction method consists of a change in the transport algorithm that reduces the statistical error (i.e. the variance) in the dose distribution in a way that is unbiased, i.e. that does not change the average value of dose in each voxel. For a prescribed level of accuracy xcex5, variance reduction has the effect of lowering the computational time Tcomp. Conversely, a variance reduction method allows a better (i.e. smaller) accuracy level xcex5, for the same computational time Tcomp.
Some of these variance reduction methods work by splitting particles into multiple representative particles at certain points in the transport simulation. The increase in the number of particles is balanced by assignment of a weight to each of the particles. The reason for employing particle splitting and weights is to optimize the use of the particles. For simple interactions (i.e. cross sections with little variation), relatively few particles are needed to get an accurate sample of the effect of the interaction. For complex interactions (i.e. cross sections with a lot of variation) many particles are needed to adequately sample the interaction. So the simulation is performed most efficiently if few particles are used for the simple interactions, but the few particles are split into many when a complicated interaction occurs.
In spite of its accuracy, EGS4 has not been widely implemented in clinical applications, because it may take as long as 20 or more hours for a single dose computation on a modern workstation, even with use of the variance reduction methods in [7]. In a typical cancer center, however, the dosimetrist who is responsible for creating the treatment plan can only spend about an hour or so with each patient, and creation of the treatment plan requires at least several dose computations. Embodiments of the present invention, NXEGS, remove this obstacle by acceleration of the Monte Carlo code to clinically acceptable speeds.
In addition to simulation of particle transport and dose deposition within a patient or phantom, the radiotherapy treatment planning also requires determination of the radiation emitted from the linear accelerator and its modification by the beam modifiers. The beam modifiers, shown in FIG. 1.2, are adjacent to the linear accelerator""s exit plane (located at the exit of the linac) and modify the radiation beam emitted from the linear accelerator.
BEAM
The program BEAM, developed at NRCC [3], is an EGS4 user""s code for simulation of linac output and its modification by the beam modifiers. It uses methods similar to those in EGS4, and for a number of linear accelerators, BEAM has been shown to reproduce the beam output quite faithfully (e.g. see [4]). On the other hand, BEAM requires complete and accurate knowledge of the components of the linear accelerator (e.g. geometry and constitutive materials), which is not readily available from the manufacturer. In contrast, embodiments of the present invention, NXEGS, includes a beam commissioning tool that uses dosimetry calibration data, rather than manufacturer specifications, to construct a source model that faithfully represents the output of the linac.
BEAM also simulates the radiation transport through the beam modifiers. The current implementation is incomplete, however, in that it does not include dynamic treatment modifiers, such as dynamic wedges and multi-leaf collimators (MLC), and it does not include compensating filters. In addition, simulation through the beam modifiers in BEAM is not directly connected to a simulation method for the patient and phantom. Also, BEAM may take many hours to simulate transport through a set of realistic beam modifiers. Embodiments of the present invention, as implemented in the computer code NXEGS, include the full range of available beam modifiers, including dynamic wedges, MLC and compensating filters. Moreover, NXEGS is directly connected to transport simulation and dose computation in the patient/phantom, and it performs computations of transport through the beam modifiers in a clinically acceptable time.
Other Programs
Other related methods include MCNP and Peregrine. Similar to EGS4, MCNP is a general method for Monte Carlo simulation of particle transport. It has been applied to radiation therapy, e.g. [5], with results that are similar to those of EGS4, but its use is not as widespread as that of EGS4.
Peregrine refers to a project at the Lawrence Livermore National Laboratory (LLNL) on Monte Carlo simulation for radiotherapy [6]. It includes a wider range of particle types [14] than EGS4, and uses updated collision cross sections. Its beam commissioning approach [13] is similar to that of BEAM, in that it performs direct simulation of the linac and requires full specification of the internal components of a linac. Peregrine applies parallel computation to accelerate the simulation of particle transport in the beam modifiers and patient.
Several other techniques have been developed for acceleration of Monte Carlo simulations for particle transport. These methods reduce the statistical error (variance) in the dose distribution produced from simulation of N source particles in the radiation beam. This reduces the value of N that is required to meet a given accuracy goal, and so it accelerates the computation.
The Voxel Monte Carlo (VMC) method of [17] and the Dose Planning Method (DPM) of [18] completely change the basic particle transport method in EGS4 in order to accelerate the computation. The VMC method uses precomputed transport steps to increase the size of the steps. The DPM method reorders the combination of particle motion and particle interactions.
Linear digital filtering techniques have been proposed in [15] and implemented in [16] for reducing the statistical error in the dose distribution after completion of the transport simulation. These filtering techniques remove any high wavenumber components in the dose distribution. Earlier application of filtering to inverse dose calculation for radiation therapy appeared in [19].
The new acceleration methods of embodiments of the present invention are distinct from these earlier methods. Embodiments of the present invention use new variance reduction methods and a new nonlinear image processing method, which are unrelated to the parallel computing, modified transport and linear filtering methods described above.
An embodiment of the invention is implemented as two software tools: a beam commissioning tool and a Monte Carlo dose calculation tool collectively referred to as NXEGS.
Beam commissioning is the process of obtaining a radiation source model to represent a specific linear accelerator in the clinic. The radiation source model is then used as input for the Monte Carlo dose calculation tool.
Linear accelerators can be run as either electron beams or photon beams; the calibration is performed separately for these two modalities. The energy (MeV) of the beam must also be specified.
For each beam type (electron or photon; both treated as particles) a large data set is obtained by measurement of dose deposition in a target (the xe2x80x9cphantomxe2x80x9d) consisting usually of water. The required data set consists of scanning data and non-scanning data for a plurality of values of beam energy and beam width. Scanning data consists of dose values along the central axis of the beam or along a line in a cross section of the beam. Non-scanning data consists of dose values at specified points on the central axis of the radiation beam. Each of these is performed for beams of various widths and energies. Dose deposition values are typically measured in the phantom with a small ion chamber. This data collection requirement is similar to that required for existing commercial RTP systems, which are non-Monte Carlo implementations, such as FOCUS, sold by Computerized Medical Systems, St. Louis, Mo.
From the measured data, the commissioning process produces a source model that accurately represents the output of the actual linac. The source usually represent the phase space in terms of particles originating from several independent subsources; e.g. a focal point source at the beam center, extra-focal sources, and contamination sources. The radiation beam from the linear accelerator is simulated in the Monte Carlo transport computation by sampling particles from the subsources. For each such particle, the initial position, direction and energy of the particle is chosen at random by sampling from the subsources, which can be described as a phase space distribution (PSD) for the particles. The PSD consists of the probabilities for each combination of values of the initial positions, direction and energy for emitted particles. The initial position may be specified as a point on a sampling plane. For example, it may be taken to be at a source point within the linear accelerator, at a point on the exit plane of the linear accelerator, or in a plane at the top of the scrapers. This PSD can be sampled to obtain an arbitrary number of distinct particles which can then be used in an RTP simulation.
The first step in the dose calculation is to setup the geometry and material properties of the phantom or patient and beam modifiers. The patient geometry includes both the patient and any bolus (a treatment aid placed on the patient to minimize the effects of body curvature). A geometric grid is setup with a prescribed spacing (e.g., 5 mm). The outer boundary of the grid is its xe2x80x9cscope.xe2x80x9d Once a particle leaves scope, it is no longer followed. Material properties of the patient are input as CT scan values, which are converted into material identification. The properties (e.g., interaction cross sections) of the identified material are, in turn, described through the PEGS4 (Preprocessor for EGS4) data set [1]. The PEGS4 data set consists of cross-section data that describes the statistics of the outcome for each of the interactions included in the EGS4 simulation method, for each energy of the incident particles, for each of the materials found in radiotherapy, and for each density of those materials.
The beam modifiers are used to shape the intensity profile of the radiation beam that is emitted from the linac. These include wedges, cutouts, ports, blocks, multi-leaf collimators (MLC) and compensating filters. For a given source model, particles are emitted from the linac exit plane or some other sampling plane, then transported through the beam modifiers. Transport through MLCs and custom ports is handled through modeling; transport through the other beam modifiers is simulated by Monte Carlo. Exiting particles from the beam modifiers are collected on a xe2x80x9ctrapxe2x80x9d from which they are then emitted into the patient geometry. The geometry of the phantom and that of the beam modifiers have different centers and coordinate axes.
The basic transport simulation and dose deposition algorithm in NXEGS is the same as that of EGS4. However, NXEGS includes a new implementation of the EGS4 algorithms as a C++, object-oriented software system and with new features and improvements described herein.
Acceleration of the RTP simulation to clinically acceptable speeds is achieved in NXEGS through new transport methods, variance reduction and postprocessing methods. These methods reduce the statistical noise for a fixed number of source particles; equivalently they reduce the number of particles required for a prescribed accuracy level. Since computational time remains roughly proportional to the number of particles, this speeds up the computational time.
Variance reduction methods change the transport algorithm in order to reduce the generation of statistical noise. EGS4 employs a number of variance reduction techniques, including range rejection, interaction forcing, particle splitting and Russian roulette [7]. NXEGS uses a new particle splitting technique for the electrons affected in Compton scattering. One more useful addition is the distributed deposition from the xe2x80x9cdeadxe2x80x9d electrons (those electrons whose energy falls below the cutoff), which we refer to as xe2x80x9czombiexe2x80x9d transport.
Postprocessing in NXEGS is applied after the transport simulation is completed to reduce the statistical noise in the resulting dose distribution, without changing the algorithm for particle transport. This postprocessing is a form of image processing, in which the dose distribution values define the image. It is a new method for Monte Carlo acceleration. It is nonlinear, in that it first extracts information on the statistical error; then uses that information to eliminate some of this error from the distribution.
In summary, NXEGS provides a complete solution to beam commissioning and dose computation for RTP applications. Using standard dosimetry data, it produces a source model to match the radiation output from a given linac. For the beam modifiers, NXEGS handles the full range of beam modifiers, including dynamic beam modifiers, MLC and compensating filters. While the basic transport and dose deposition method is effectively identical to that of EGS4, novel acceleration methods in NXEGS greatly speedup the Monte Carlo simulation to achieve clinically acceptable computation times with high accuracy in the dose distribution.
Embodiments of the invention include a method of commissioning a radiation source which inputs measured dose data into a data processor. The measured dose data is derived from exposing a phantom to radiation from the source; and measuring the radiation dose to obtain a measured dose in the phantom resulting from the exposing step. The measured dose is measured at a plurality of points within the phantom. At least some of these points are axial points located at positions along a substantially central axis of the radiation source and others of these points are transverse points located at positions along an axis transverse to the central axis. After inputting he measured dose data, one performs a Monte Carlo simulation of the radiation source to determine a simulated dose at the plurality of points; and then models the radiation source using the simulated dose and the measured dose.
The dose data is obtained from positioning a phantom to receive radiation from the source, exposing the phantom to radiation from the source, and measuring the radiation dose in the phantom resulting from the exposing step. The number of dose measurement points may be 10 or greater and is preferably 50 and most preferably is a large number of dose measurements on the order of 100 points. The above described method has direct application to for use in radiation therapy planning.
According to other embodiments of the invention, there is described a method of modeling a radiation source which emits radiation in the form of particles. One step of the method includes defining at least one simulated sub-source comprising one of a focal source and an extra focal source; each simulated sub-source having a relatively small interval of energy values and a relatively small interval of angular values. Further, one performs a Monte Carlo simulation of particles leaving the at least one simulated sub-source to determine a simulated dose distribution; and then determines a phase space distribution of radiation leaving the at least one simulated sub-source using said simulated dose distribution and actual dose measurements.
According to other embodiments of the invention, there is described a method of performing radiation planning therapy in which one performs a CT scan of a patient to provide patient-dependent information in a region of the patient to be treated; utilizes a source model based on dosimetry data which includes at least 10 dosimetry points, and typically about 100 points, provides treatment head information concerning the characteristics and geometry of a treatment head; and simulates a dose distribution using a Monte Carlo calculation based on a source model, the patient-dependent information and the treatment head information.
Embodiments of the invention are characterized by a method of defining a source model of an accelerator emitting radiation. This method measures the actual radiation dose distribution from the accelerator over a spatial volume defined within a phantom and for a plurality of accelerator energies; uses a Monte Carlo simulation to calculate a transport matrix which relates the actual radiation dose distribution to a source phase space distribution, where the source phase space distribution is defined as the probability distribution of the position, energy and direction of radiation from the accelerator. Finally, the method calculates the phase space distribution using the transport matrix and the actual radiation dose distribution.
Yet other embodiments of the invention use a method of simulating radiation transport through a treatment head of a linear accelerator in which the treatment head has at least one beam modifier. The method employs inputting into a data processor, parameters corresponding to a physical description of the treatment head and the at least one beam modifier included in the treatment head; simulating the introduction of particles from a source into the treatment head; and performing simulated transport of particles through the at least one beam modifier in the treatment head through at least one of a Monte Carlo method or a throughput function.
There is further described according to other embodiments of the invention, a method of radiation treatment planning which includes performing one of (1) commissioning a linear accelerator to obtain a source model or (2) using a source model from a linear accelerator that has already been commissioned, wherein the commissioning includes utilizing at least 10 dosimetry points; simulating one or more beams based on a source model; specifying a simulated treatment head, including (1) geometry and (2) material properties of at least one beam modifier for each beam; setting up a simulated patient or a simulated phantom; simulating particle transport through the simulated treatment head using at least one of a Monte Carlo method or a throughput function for each beam; simulating particle transport through the simulated patient or simulated phantom using the Monte Carlo method; calculating a radiation dose in the simulated patient or simulated phantom; and providing an output of the calculated dose. In addition to specifying a simulated treatment head, including (1) geometry and (2) material properties, one may also specify motion of the at least one beam modifier.
Embodiments of the invention are further described as a method of specifying at least one beam modifier for use in Monte Carlo simulation of a particle motion through a treatment head which includes a radiation source and at least one beam modifier in which the method includes inputting into a data processor a geometry specifying said at least one beam modifier by performing at least one of the following three steps: (i) for a beam modifier in the shape of a polyhedron, inputting the coordinates of each vertex of the polyhedron and specifying the edges of connecting vertices, and inputting the position and orientation of the polyhedron; (ii) for a beam modifier in the form of a multi-leaf collimator, inputting the number and size of each leaf; and (iii) for a beam modifier in the shape of a compensating filter, consisting of a flat side and a curved shape of arbitrary complexity opposite the flat side, inputting the height of the curved shape above each point of the flat side, and specifying whether the flat side points toward or away from the radiation source. Finally, the method includes inputting a type and density of material that constitutes the at least one beam modifier.
The above method may include movement of the at least one beam modifier. In such a case, one further specifies movement of the at least one beam modifier by inputting a series of frames defined by specifying at least one of (1) positions and (2) orientations of the at least one beam modifier and a particle weight corresponding to each of the at least one of the frames; randomly selecting a frame for each of a plurality of simulated particles emitted from the radiation source; and applying the corresponding particle weight to each of the plurality of simulated emitted particle for use in the Monte Carlo simulation. In this manner, utilizing the particle weight simulates a time interval for each of the frames, and thus simulates the change of the at least one of the position and orientation of the at least one beam modifier so as to simulate moving the at least one beam modifier during the Monte Carlo simulation.
According to other embodiments of the invention, there is described a method for computing the dose deposition due to particles whose energy is below a cutoff, as part of a Monte Carlo simulation. The method comprising the steps of determining that a given particle has reached a cutoff energy within a given voxel; and for particles below the cutoff energy, performing the steps of: extending the particle path in a straight path from the point at which the cutoff was reached in the direction of the particle motion at said point; and simulating a constant energy deposition rate using a rate constant along the straight path .
Other embodiments of the invention are directed toward a method of simulating the transport of electrons using a Monte Carlo calculation. For a given energy of a simulated particle, one determines a step size in the Monte Carlo calculation to correspond to a fixed number of elastic scatterings wherein said step size is between an initial and end position and the particle path is a straight path between said initial and end positions. One further determines whether a discrete event occurs along said particle path based on the cross section of the particle within the material through which the particle path traverses. If a discrete event does not occur along the particle path, then the method utilizes the steps of: randomly selecting a first intermediate point between the initial position and the end position along said particle path; determining a scattering angle at the first intermediate point based on a tabulated cross section of the material existing at the first intermediate point; and simulating a Molixc3xa8re scattering event at the first intermediate point with the determined scattering angle wherein the Molixc3xa8re scattering corresponds to the fixed number of elastic scatterings. The scattering angle determines a new end point for the particle path.
Still further embodiments of the invention are directed toward a method for commissioning a source model for an electron beam wherein the source model includes a focal electron subsource and a focal photon subsource, both lying on the central axis of a linear accelerator, wherein said source model includes a sampling plane perpendicular to the central axis in which simulated electrons and photons are specified by their position r from the central axis and their energy E, wherein the position r defines an angle theta, and the position r is randomly selected and the particle is assigned a weight through fluence functions for the electron and photon subsources and the energy E is determined through energy spectra for the electron and photon subsources. The method comprising the steps of: selecting an initial estimate for (1) positions of the electron and photon subsource along the central axis, (2) parameters defining the energy spectra for the electrons and photons; and (3) parameters defining the fluence functions for the electrons and photons. One then uses the initial estimate and simulates a plurality of dose distributions, including output factors, percentage depth dose (PDD) curves, and dose profiles. One then compares the simulated output factors to the measured output factors and determines an improved estimate for the positions of the electron and photon subsource along the central axis; compares the simulated PDD curves to the measured PDD curves and determines an improved estimate for the parameters defining the energy spectra for the electrons and photons; and compares the simulated dose profiles to the measured dose profiles and determines an improved estimate for the parameters defining the fluence functions for the electrons and photons.
Another embodiment of the invention concerns a method for defining a source model for an electron beam wherein the source model includes the following: (a) defining a focal electron subsource lying on the central axis of a linear accelerator; (b) defining focal photon subsource lying on the central axis of the linear accelerator; (c) defining a sampling plane perpendicular to the central axis ; (d) specifying fluence functions, which are functions of the position r from the central axis, for each of electrons and photons; (e) specifying energy spectra for each of the electrons and photons; (f) determining the direction of each of the electrons and photons by the straight line direction from the respective focal subsource location to a point in the sampling plane; and (g) specifying a scattering of electrons by varying the direction of the electrons from the determined direction from step f).
Another embodiment of the invention is directed toward a method for commissioning a source model for a photon beam wherein the source model includes a focal photon subsource and a focal electron subsource, both lying on the central axis of a linear accelerator, and an extrafocal photon subsource lying on a plane perpendicular to the central axis, wherein said source model includes a sampling plane perpendicular to the central axis in which simulated electrons and photons are specified by their position r from the central axis and their energy E, wherein the position r defines an angle theta, and the position r is randomly selected and the particle is assigned a weight through fluence functions for the electron and photon subsources and the energy E is determined through energy spectra for the electron and photon subsources, The method comprising the steps of: (a) determining the spatial distribution of said extrafocal photon subsource to fit output factors in air; (b) performing Monte Carlo simulation of radiation transport for a plurality of beamlets, each of which corresponds to subset of energies and angles for said focal photon subsource or the focal electron subsource, and obtaining simulated dose values at each point of a phantom for which there are measured dose values; and (c) using a matrix of said simulated dose values in determining the spectrum and fluence for the focal photon subsource or the focal electron subsource to fit the measured dose values.
According to further embodiments of the invention there is disclosed a method for defining a source model for a photon beam wherein the source model includes the following: (a) defining a focal photon subsource lying on the central axis of a linear accelerator; (b) defining a focal electron subsource lying on the central axis of the linear accelerator; (c) defining an extrafocal photon subsource lying on a plane perpendicular to the central axis of the linear accelerator; (d) defining a sampling plane perpendicular to the central axis; (e)specifying fluence functions, which are functions of the position r from the central axis, for each of the electron and photon subsources; (f) specifying energy spectra for each of the electron and photon subsources; (g) selecting a point in the sampling plane for each electron from the focal electron subsource, each photon from the focal photon subsource and each photon from the extrafocal photon subsource; (h) selecting a point on the plane of the extrafocal photon subsource; (i) determining the direction of the photon from the focal photon subsource and the electron from the focal electron subsource by the straight line direction from the respective focal subsource location to the respective point in the sampling plane; (j) determining the direction of the photon from the extrafocal photon subsource by the straight line direction from the point on the plane of the extrafocal photon subsource to the respective point in the sampling plane; and (k) specifying a scattering of electrons by varying the direction of the electrons from the determined direction from step i).
Yet another embodiment of the invention is directed toward a method of error reduction of Monte Carlo simulation radiation transport. This method comprising the steps of: (a) calculating a dose distribution fS(x) using a first Monte Carlo simulation using a number, n, of particles and a random number sequence; (b) calculating a dose distribution gS(x) based on a simplified representation of the radiation transport using a second Monte Carlo simulation with said number of particles n and said random sequence; (c) calculating a dose distribution gL(x) based on said simplified representation using a third Monte Carlo simulation with a number of particles m where m is greater than n; and (d) reducing the error in the dose distribution fS(x) by using fS(x)xe2x88x92gS(x)+gL(x) as the error reduced dose distribution.
More generally, the above described embodiment concerns a method of error reduction in radiation transport which performs a Monte Carlo simulation transport calculation to determine a dose distribution in a region of interest; and reduces the error of the dose distribution using control variates.
According to another embodiment of the invention, there is taught a method of accelerating computation of dose distributions in a target area using the steps of: (a) mathematically constructing a fine grid throughout regions of the target area where the density of the target area is non-uniform; (b) mathematically constructing a coarse grid throughout regions of the target area where the density of the target area is uniform; and (c) performing a Monte Carlo simulation of the dose distribution using transport steps defined within the coarse and fine grids wherein in the coarse grid, larger transport steps are used relative to transports steps used in the fine grid.
Yet another embodiment of the invention involves a method of improving the accuracy of a Monte Carlo simulation in radiation treatment planning for particle transport within a target region of interest. This method comprising the steps of: (a) selecting an end point position within the target region of interest; (b) selecting an initial position outside of the target region; (c) selecting at least one intermediate point on a particle path connecting the initial and end positions randomly chosen using a bi-directional probability distribution; (d) repeating step c) for multiple particle paths having different intermediate points; and (e) calculation the total dose within the target region resulting from the summation or average of the doses of the multiple paths.
Embodiments of the invention are further described as a method of improving the simulation of a Compton scattering event in computing dose distributions of radiation resulting from Compton scattering events. This method comprising performing a Monte Carlo transport simulation using an electromagnetic source; determining Compton scattering events along a simulated transport path of the radiation in a treatment region of interest; and at each Compton scattering event, simulating the generation of a plurality of electrons distributed along the transport path.
Another embodiment of the invention includes a method for error reduction of Monte Carlo simulation for radiation transport by postprocessing the resulting energy or dose distribution using at least one of nonlinear filtering or image processing techniques. The radiation transport may be applied to radiation therapy and the may utilize one of the following techniques: (a) a Donoho-Johnstone soft thresholding; (b) a Osher-Rudin equation; and (c) a local cosine method.